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ABSTRACT 

We present a numerical implementation of the renormalization group (RG) for 
partial differential equations, constructing similarity solutions and travelling 
waves. We show that for a large class of well-localized initial conditions, succes- 
sive iterations of an appropriately defined discrete RG transformation in space 
and time will drive the system towards a fixed point. This corresponds to a 
scale-invariant solution, such as a similarity or travelling-wave solution, which 
governs the long-time asymptotic behavior. We demonstrate that the numerical 
RG method is computationally very efficient. 
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1. Introduction 

Recently, renormalization group (RG) theory has been apphed to a number of 
non-equihbrium physical systems"^ without any statistical aspect, to study the 
long-time or large-scale intermediate asymptotic behavior. In particular, the im- 
portant cases of similarity solutions"^ ^ u(x^t) = t~"f (xt~^^^ or travelling- wave 
solutions"^*^ u(x^t) = f (x — vt) were treated. In all of these cases, generically 
referred to as exhibiting intermediate asymptotics of the second kind,"^ the expo- 
nents a, fi or the velocity v cannot be determined by simple dimensional analysis 
or from conservation laws, but are determined only by solving the full prob- 
lem itself. Previous work demonstrated that these exponents are the anomalous 
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dimensions of the field theoretic RG, and may be calculated using RG 

to remove systematically divergent or secular terms from a naive perturbation 

expansion. 

Although RG is often formulated in connection with perturbation theory, it is 

essentially non-perturbative and has a direct geometrical interpretation, as Wilson 
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showed in the context of critical phenomena. The purpose of this present paper 
is to explore the long-time asymptotic behavior of certain simple non-equilibrium 
systems within this geometrical picture. We will see that the fixed points of the 
appropriately defined RG transformation correspond to scale-invariant solutions, 
such as similarity or travelling- wave solutions. The simple numerical method 
given here, which exploits this picture, extends the practical applicability of the 
RG to problems where there is no small perturbation parameter, in contrast to 
our previous perturbative work. 

Other numerical methods, based on rescaling and thus closely related to the 
RG, have been applied to study (e.g.) finite time blow-up in the two-dimensional 
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nonlinear Schrodinger equation and in axisymmetric three dimensional Eu- 
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ler dynamics. Numerical RG methods, particularly Monte Carlo RG, are widely 
used in equilibrium statistical mechanics to calculate equilibrium phase diagrams 
and anomalous dimensions.^ Monte Carlo RG has also been applied to a partic- 
ular non-equilibrium system — the driven diffusive gas — to calculate scaling 
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exponents of the nonequilibrium steady states. 

We begin by considering numerical RG transformations to obtain self- 
similarity solutions. A RG transformation i?;,,/?, which depends on two parame- 
ters, a dilation parameter b and an exponent /3, is defined on the space of functions 
u(^ x ^ at some arbitrary time t by u'(x^t) = Ri,f]{u(x^t)}. This transformation 
involves the following three steps: 

(1) Evolve the function u(x^ t) forward over finite time, from t to t' = 6t, 6 > 1, 
using the governing PDE, and call the result u{x, t'). This step is analogous 
to the block spin transformation used in critical phenomena. 

(2) Rescale x by defining x' = b~^x^ so that u(x',t') = u{x'U\t'). 

(3) Rescale the function u itself by an amount Z{h) = u(0^t')/u(0^t)^ so that 
u'{x',t') = Z{b)u{x',t'). 

Thus, the resulting RG transformation yields u'(x^t) = Z (b)u (^b^ x ^ bt^ . The 
basic idea is that any fixed point u* = R{u*} corresponds to a similarity so- 
lution. For the simple systems discussed in this paper, only one fixed point 
exists, which is a stable attractor. We will see that repeated iterations of the RG 
transformation drive the function u to its fixed point value more rapidly than 
straightforward time evolution for large times. The semi-group property of this 
RG transformation (6 > 1) : Ri,^ f]Ri,^ f] = Ri,^i,^ f]^ simply implies Z(b) = b" or 
equivalently a = dlog Z / dlogb. The parameter /3 must be varied to search for 
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fixed points of this RG transformation, which may not exist for arbitrary /3. For 
appropriate values of /3, after an infinite number of successive iterations of the 
RG transformation defined above, the fixed point, wiU be attained, if 

the initial conditions are in the basin of attraction of the fixed point, i.e., we 
have u*(x^t) = h^"u (fc^^'x, , as n — )• oo. By setting 6" = A/t, we find for 
any constant A, u*{x^t) = t~"(f> (^xt~^^^ as t — )• oo, where ^ is a scaling function 
to be determined. In principle, the RG transformations outlined above - either 
continuous or discrete - can be implemented in many different ways, but here we 
will use a discrete formulation, which is more useful for numerical applications. 

2. Numerical Method for Similarity Solutions 

In our numerical implementation, we adopt a finite-difference Euler scheme 
to discretize the continuous version of the PDE at hand. We write t = zAt, x = 
j Ax, z = 0, 1, 2, • • • , M, j = 0, 1, 2, • • • , iV, with At a time step and Ax a mesh size. 
By simply evolving the PDE over a short, finite time from t = to to ti = bto with 
no = (6— l)to / At time steps, we obtain a set of data \^u'{x = jAx, ti = 6to)} . The 
step (2) is best realized by rescaling the mesh size Ax, rather than by changing 
the discrete sites j = 0, 1,2, • • • ,iV at each iteration. Thus, after one iteration, 
the new mesh size is {Ax)i = Ax. Finally, we implement step (3) by rescaling 
the whole discrete set of data by a factor Zi{h) = u'(0^bto)/u(0^to). Thus, after 
the first RG iteration, a set of coarse-grained data is obtained at fixed sites j: 

uW(xi,to) = Zi{b)u'{b^xi,bto), (2.1) 

with xi = j(Ax)i^j = 0,1,2,---. This data set will constitute the new initial 
condition for the second RG iteration. The time step At and the total number 
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of time steps no = (b — l)to/At are kept fixed through all the RG iterations, so 
that after n RG iterations, we will have calculated the n-th iteration 

to) = Zn(b) ■ ■ ■ Zi(b)y' 6"to) , (2.2) 

where Xn = j{Ax)n with a rescaled mesh size (Ax)n = b~"'^Ax. Thus, in our 
numerical RG transformation, the equation is always solved over a finite time 
interval, and the system size N(Ax)n continues to shrink with the number of RG 
iterations increasing. This behavior can be understood by noting that the limit 
x/t^ ^ in a possible similarity solution u(x^t) = t~" (f>(x /t^) can be realized 
by letting x (or mesh size Ax) — )• and keeping t fixed at a finite value, instead 
of letting t — )• oo and keeping x fixed. Therefore, the numerical RG method 
is conceptually different from direct numerical integration (DNI), in which we 
simply evolve the PDE over a very long time t and generally have to choose a 
very large-size system. 

In practice, a potential difficulty is the choice of parameter /3. For the fixed 
point, /3 is determined by requiring that, if the curve t"u{x^ t) vs. xt~l^^ is plotted, 
all the data should collapse onto a single curve, namely the universal scaling 
function cf). A simple alternative is to require that, if \ogu{x = 0,t) vs. logt is 
plotted, all the data should lie on a straight line, enabling one to read off the 
slope —a. In practice, the correct value fi is selected in the following way. Given 
an estimated interval for /3, /3i < /3 < /32, we calculate in the last RG iteration, 
the slopes d\ogu{x = O,t)/c/logt at t = to, ai, and at t = 6to, a2, respectively. 
For arbitrary values of /3, the difference between ai and a2, should be a function 
of /3, i.e. Aa = a2 — ai = f{/3). Only for certain /3, will Aa be equal to zero. 
This is a typical root finding problem, if we regard /3 as the root of the function 
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/. The appropriate value of /3 = /3* is then easily found using any given root 
finding algorithm. We have found this method to converge rapidly, and still be 
faster than DNI. 

The first example we consider is the so-called modified porous-medium equa- 
tion 

dtu = DAu^+'', (2.3) 

where A represents a (/-dimensional Laplacian, n is arbitrary, D = 1 for dfU > 

5 

and D = 1 -\- e for dfU < 0. This problem has previously been treated 
analytically using the RG method, for e ^ 1. We have applied the numerical 
discrete RG transformation to several special cases, including the problem of 
gravity-driven groundwater (n = l,c/ = 2). For simplicity, we will take the 
Barenblatt equation (n = 0,c/ = 1) to illustrate our basic ideas on numerical 
RG. This equation describes the pressure u during filtration of an elastic fluid in 
an elasto-plastic porous medium. For e 7^ 0, the material exhibits hysteresis and 
the diffusion decays anomalously, with the long-time asymptotics of the similarity 
form u{x, t) ~ t~^^^'^~^"^ f{xt~^^'^ , e), where the anomalous dimension is calculated 
using perturbative RG^ to be a = e/(27re)"^/^ -|- O(e^). 
A simple, explicit^ discrete scheme can be written as 

z + 1) = z) + DrAuij, z), (2.4) 

where u{j,i) = u{x = jAx,t = zAt),r = At/(Aa;)^, the discrete Laplacian is 

Au(j, z) = u{j + 1, z) - 2u(j, z) + u{j - 1, z), (2.5) 

and D = 1 if Au(j^ z)>0,-D = l-|-eif Au(j^ z) < 0. The condition for numerical 
stability Dr < 1/2 must be imposed in this explicit scheme. The numerical values 
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of the parameters 6, At, Ax are chosen in such a way as to not only resuh in the 
quickest rate of convergence to the true solution, but also to attain sufficient 
accuracy. Typical values of At and Ax are chosen as 0.01 and 2.0, respectively, 
because a too small value for r = At/{AxY gives a high accuracy, but needs 
too large a number of time steps no, whilst a too large value for r (which still 
satisfies the numerical stability condition) is not accurate enough. Although 
in principle any value b > 1 can be chosen, we choose a dilation parameter b 
typically about 1.02. For too large a value of 6, large number of time steps no 
would be required, and the system would shrink too quickly even after only a few 
number of RG iterations (the mesh size at n-th round is (Ax)n = Ax/b""^). For 
a very small 6, a large number of RG iterations are necessary to produce good 
enough accuracy, and this takes too much computer time. Thus, we must adopt 
a compromise between the number of RG iterations and the size of the finite time 
interval used in each iteration. The typical number N of RG iterations is chosen 
as 100 < iV < 500, depending upon how accurate the solution is required to be, 
whilst satisfying the condition for numerical stability DAt/(Aa;)^ < 1/2 even 
after N RG iterations. 

In figure 1, we show how the anomalous exponent a extracted from Z(b) varies 
with the number of RG iterations N and converges to the value a = 0.6975 in 
the case e = 1.0, starting from different localized initial conditions, where we 
have taken only a small number of time steps no. As we see, an advantage of 
our numerical RG is that, even if we choose a very small number of time steps 
no {e.g. no = 10) in each RG iteration (an obviously poor approximation), after 
a sufficient large number of RG iterations, the final result for the exponent a is 
still very accurate. 
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For purpose of comparison, we also perform a conventional DNI of the Baren- 
blatt equation in d = 1 and 2 dimensions using the same explicit scheme. We 
choose the same initial conditions as in the numerical RG calculation, Ax = 2.0 
and At as large as possible, whilst satisfying the numerical stability condition 
r = At/{ AxY < l/2d. We compare both CPU times consumed on the same com- 
puter and the total number of iterations performed so that the same accuracy of 
the solution {e.g. the same relative error) is reached. 

Suppose that the number of time steps taken to have the same accuracy 
is ni and n2, respectively, in the DNI and the numerical RG method, and the 
number of RG iterations is N . An initial condition with compact support will 
be non-zero at most over n -\- 1 additional spatial grid points after the n-th time 
iteration. Therefore, the total number of iterations is estimated to be of order 
Ni = 0{n'j^~^^) and N2 = 0{Nn2~^^) for ni,n2 ^ 1, respectively, in d dimensions. 
If ni/n2 ^ 1, then Ni ^ iV2, which implies that DNI is much slower than 
the numerical RG method. To achieve this goal, we choose N typically about 
100 — 500, and n2, typically about 20 — 50. 

Depending upon the initial conditions, the explicit numerical RG scheme is 
about 5 — 10 times faster than DNI in d = 1 dimension, and about 50 — 100 
times in d = 2 dimensions, where the requisite accuracy is 1 — 2%. We find 
that, in DNI, the more unsymmetrical and unsmooth the initial conditions are, 
the more slowly the solutions converge, especially when a very high accuracy 
is required. It is not unreasonable to expect that our numerical RG method 
will be about several hundred times faster than DNI in d = 3 dimensions, with 
arbitrary, unsymmetrical, localized initial conditions. Therefore, it seems that 
the numerical RG method is computationally more efficient than DNI of equation 
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of motion, because the RG transformation effectively coarse-grains uneven parts 
of the initial data, and drives more quickly all the initial conditions in its basin 
of attraction towards its fixed point. Another advantage of our numerical RG 
method is that, because only a small number of time steps and a small size system 
are needed, there is no difficulty at all in applying it to two or three dimensional 
systems, whilst if DNI is used, in order to obtain a reasonably accurate solution, 
a huge number of time steps and a large system is required. 

We have also applied this numerical RG to study the universal behavior of 
the long-time asymptotics. We choose a number of different initial conditions, 
sufficiently localized and integrable, i.e., with J dx u(x, 0) = Qo bounded, such as 
a Lorentzian, step function, and Gaussian, or even those initial distributions with 
a major peak at origin and other minor multiple peaks elsewhere. In all cases, as 
long as the initial conditions are well-localised, the numerical RG yields the same 
anomalous exponents and scaling functions for the long-time asymptotic dynam- 
ics, although the rates at which they converge are different. It is also found that 
the long-time asymptotic behavior is independent of specific numerical schemes, 
explicit or implicit; nevertheless, the explicit ones are less time-consuming than 
the implicit ones, and for this reason we prefer to adopt the explicit scheme in 
all our simulations. 
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3. Numerical Method for Travelling Waves 



In this section, we present a numerical RG transformation for travelling waves 
interpolating between stable and unstable states. In earlier work, we applied a 
variational principle and perturbative RG to study the dynamical velocity selec- 
tion mechanism,"^*^ where it was proposed that physically relevant travelling-wave 
solutions must be structurally stable in a precise sense. 

We define a RG transformation i?j „ on the space of functions u(x^t) at 
some time t by u'(x,t) = Ri,^j;{u(x,t)}, depending on two parameters, a dilation 
parameter b and a speed v. Two steps follow in this transformation: 

(1) Evolve the function u(x^t) forward over a finite time, from t to t' = (b -\- 

6 > 0, using the governing PDE, and call the result u'(x^t'). 

(2) Rescale x by defining x' = x -\- vbt, i.e. shifting x by an amount of dis- 
placement v{t' — t) = vbt, where v must be varied to search for the fixed 
point. 

Unlike the similarity case, there is no need to rescale the function u itself, 
because the rescaling factor Z{h) is actually equal to one (This corresponds to 
the exponent a = if we transform the travelling-wave form f{x — vt) into the 
similarity form T~" f(XT~^) by T = logt,X = logx.) Thus, we have 

v'(x,t) = u(x + vbt,(b + l)t). (3.1) 

This transformation is useful in that any fixed point of i?j „ is a travelling- wave 
solution. This RG transformation generates a semi-group (6 > 0): Ri,^ vRh2 v = 
Rl,^-\-b2 V A. fixed point will be reached, after an infinite number of RG iterations, 
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if the initial conditions are in the basin of attraction of it, i.e., 



u*(x^ t) 



u{x -\- vnht, {nh + l)t) , 



n 



oo. 



(3.2) 



By translationaUy shifting the solution by an amount v{nh + we have 



must vary the unknown velocity v to find the fixed point by requiring that 
M(0,t) = u(d^t')^ where t' = (1 + d = vo(t' — t) = v^ht, and vq is the 
desired velocity. Only when the profile is shifted back along the —x direction by 
an appropriate displacement will the shifted solution coincide with the original 
one. If we were to choose a velocity v larger (or smaller) than uq, then the solu- 
tion would be overshifted (or undershifted) back and u(vbt^t') would be smaller 
(or greater) than M(0,t). We emphasize that this is completely different from 
DNI, in that during RG iterations the unknown parameter v is being varied and 
we are taking the coarse-grained data as our new initial conditions. 

The numerical discretization procedure is almost identical with that for the 
similarity case. Due to discretization, the number of grid sites by which we shift 
back to rescale x must be an integer /, and we have to take a varying mesh 
size at each RG iteration. Ax = v/m with m a large integer number, so that 
/ = vto/Ax = mto is an integer. Typically we choose 5 < / < 20 to guarantee a 
sufficiently high accuracy of computation. Just as in the similarity case, given an 
estimated interval for < u < U2, a root-finding algorithm is used to pick out 

the dynamically selected velocity u*, by regarding the difference of u(x = 0,to) 
in the last two RG iterations as a function of Au(0^to) = f{v)^ and requiring 



u*(x^ t) 



u{x — ut), as n 



oo. As we iterate this RG transformation, we 



/(^*) = o. 
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To illustrate how this discrete numerical RG transformation scheme works, 

19 

we present Fisher's generalized population model, 

dfu = d^u + u(l — u)(l + z/m), —1<v< +00, (3-3) 

where u = 0, 1 are two steady-state points. 

When = — 1, this equation reduces to the well-known Fisher-Kolmogorov- 
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Petrovsky-Piskunov equation, for which rigorous results are known. If the 

initial condition u{x,t = 0) is assumed to be well-localized and decays as fast as 
Q-i^ for large x, then for 5 > 1, the front velocity asymptotically approaches the 
value u = 2, while for 5 < 1, the asymptotic speed is v = q-\-ll q> 2. In figure 2, 
we show numerical RG calculations for the velocity compared with the exact 
analytical result above. In doing so, we take At, Ax even as large as 0.05,1.0, 
respectively, and use only a small number of grid points to discretize t and x. 
Our simulation shows that the numerical RG transformation not only drives the 
solution towards the fixed point -travelling-wave solution more quickly than DNI, 
but also uniquely determines the correct velocity v. The typical relative error is 
only about 0.5% in this case. 

Now we move on to the general case with — 1 < < -|-oo. It is well established 
that in this case, there exists a transition, as the parameter v crosses the transi- 
tion point i^c = 2, and the corresponding velocity is v = 2 for — 1 < < 2, and 
V = {2 -\- v)/\/2v for V '>2. Starting with sufficiently localized initial conditions 
( e.g. a standard step function u{x,t = 0) = Q{—x) ), we have performed the 
numerical RG simulation, and plot our results in figure 3, in comparison with the 
exact result. It demonstrates that our numerical RG transformation does predict 
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a transition between linear-marginal-stability (pulled) and nonlinear-marginal- 
stability (pushed) cases and uniquely selects the correct front velocity for most 
natural initial conditions. 
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FIGURE CAPTIONS 



Fig. 1. The anomalous scaling exponent a as a function of the number of numerical 
RG iterations iV, in the case of e = 1.0 for the Barenblatt equation. Shown 
is the convergence to the exact value a = 0.6975 for three different ini- 
tial conditions with comparable width: Lorentzian (filled circles), Gaussian 
(filled squares) and top hat (filled triangles). 

Fig. 2. The propagation velocity u as a function of decay rate q of initial conditions. 

The points determined by the numerical RG method are denoted by •. The 
continuous curve is the exact result. 

Fig. 3. The propagation velocity v plotted as a function of v. The full curve repre- 
sents the exact result, whilst data points determined by our numerical RG 
are denoted by •. 
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